SnowUpdate Subroutine

SUBROUTINE SnowUpdate   & 
  !
  ( time, mask )       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time
TYPE (grid_integer), INTENT (IN) :: mask

!local declarations:
INTEGER (KIND = short) :: i, j, is, js
REAL (KIND = float) :: meltRate  !!snow melt rate (m/s)
REAL (KIND = float) :: refreezingRate !!refreezing rate of water retained in snowpack (m/s)
!REAL (KIND = float) :: cRefreeze = 0.05!!coefficient of refreeze (-)
!REAL (KIND = float) :: swhc = 0.1 !!snow water holding capacity (-)
REAL (KIND = float) :: alpha !!precipitation partitioning factor
REAL (KIND = float) :: liquid !!liquid fraction of precipitation (rain) (m/s)
REAL (KIND = float) :: solid !!solid fraction of precipitation (snow) (m/s)
REAL (KIND = float) :: tlower
REAL (KIND = float) :: tupper
REAL (KIND = float) :: tAir
REAL (KIND = float) :: cmelt !!snow melt coefficient (m/s/°C)
REAL (KIND = float) :: tmelt
REAL (KIND = float) :: swe_tminus1 !!swe at previous time step
REAL (KIND = float) :: ddx	!!actual cell length (m)
REAL (KIND = float) :: ds !!thickness of the saturated depth (m)
REAL (KIND = float) :: cellWidth !!cell width (m)
REAL (KIND = float) :: ijSlope  !!local slope (m/m)
REAL (KIND = float) :: qin, qout  !!input and output water discharge in snowpack
REAL (KIND = float) :: area !!cell area (m2)
!------------------------------end of declarations-----------------------------  

!check if parameter update is required
CALL SnowParameterUpdate (time)

!initialize melt grid
snowMelt = 0.


!computer inter-cell lateral water flux
QinSnow = 0.
QoutSnow = 0.
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        IF ( mask % mat (i,j) == mask % nodata ) CYCLE
        
        !set downstream cell  (is,js)
        CALL DownstreamCell ( i, j, flowDirection % mat(i,j), is, js, ddx, &
                            flowDirection ) 
        
        IF ( CheckOutlet (i, j, is, js, flowDirection) ) CYCLE
      
        !thickness of the saturated depth
        ds = waterInSnow % mat (i,j)
        
        !cell width
        cellWidth = ( CellArea (mask, i, j) ) ** 0.5
        
        !local slope
        ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx
        
        IF ( ijSlope <= 0. ) THEN
            ijSlope = 0.0001
        END IF
        
        QoutSnow % mat (i,j) = cellWidth * ds * &
                               snowConductivity % mat (i,j) *  ijSlope
      
        !output becomes input of downstream cell
        QinSnow % mat (is,js) = QinSnow % mat (is,js) + QoutSnow % mat (i,j)
        
    END DO
END DO

!loop across cells
DO i = 1, mask % idim
    DO j = 1, mask % jdim
        
        !skip nodata
        IF ( mask % mat (i,j) == mask % nodata ) THEN
            CYCLE
        END IF
        
        !potential snow melt rate
        tAir  = temperatureAir % mat (i,j)
        cmelt = meltCoefficient % mat (i,j) / 1000. / 86400.
        tmelt = meltTemperature % mat (i,j)
        SELECT CASE ( meltModel % mat (i,j) )
        CASE (1) !degree-day temperature based
            meltRate = DegreeDay ( tAir,  tmelt,  cmelt )
        CASE DEFAULT
            CALL Catch ('error', 'Snow', 'melt model not implemented: ', &
                         argument = ToString (meltModel % mat (i,j) ) )
        END SELECT
        
        ! actual snow melt, limited by available snow 
        ! in the previous time step
        meltRate = MIN ( meltRate, swe % mat (i,j) / dtSnow )
        snowMelt % mat (i,j) = meltRate * dtSnow
        
        !precipitation partitioning 
        tupper = upperTemperature % mat (i,j)
        tlower = lowerTemperature % mat (i,j)
        
        IF ( tAir <= tlower ) THEN
            alpha = 0.
        ELSE IF ( tAir > tlower .AND.  tAir < tupper ) THEN
            alpha = ( tAir - tlower ) / ( tupper - tlower )
		ELSE
			alpha  = 1.
        END IF
        
        liquid = alpha * precipitationRate % mat (i,j)
        solid  = ( 1. - alpha ) * precipitationRate % mat (i,j)
        
        ! potential refreezing rate
        refreezingRate = Refreezing ( tAir,  tmelt,  cmelt, cRefreeze % mat (i,j) )
        
         ! actual refreezing, limited by available retained water   
        ! in the previous time step
        refreezingRate = MIN ( refreezingRate, &
                               waterInSnow % mat (i,j) / dtSnow )
        
        !update snow water equivalent
        swe % mat (i,j) = swe % mat (i,j) + &
                       ( solid - meltRate + refreezingRate ) *  dtSnow
        
        !update water in snow
        area = CellArea (mask, i, j) 
        qin = QinSnow % mat (i,j) / area
        qout = QoutSnow % mat (i,j) / area
        IF ( swe % mat (i,j) > 0. ) THEN
            
            waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + &
                                      snowMelt % mat (i,j) +  &
                                      liquid * dtSnow + &
                                      ( qin - qout ) * dtSnow + &
                                     -  refreezingRate  * dtSnow 
            
            
            rainfallRate % mat (i,j) = 0.
            
        ELSE !rainfall does not contribute 
            
            waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + &
                                      snowMelt % mat (i,j) +  &
                                      ( qin - qout ) * dtSnow  + &
                                      -  refreezingRate  * dtSnow
            
             rainfallRate % mat (i,j) = liquid
             
        END IF
        
        IF ( waterInSnow % mat (i,j) < 0. ) THEN
            waterInSnow % mat (i,j) = 0.
        END IF
        
        ! check retained water does not exceed water holding capacity of snow
        IF ( waterInSnow % mat (i,j) > swhc % mat (i,j) * swe % mat (i,j) ) THEN
            !exceeding water becomes runoff
            excessWaterSnowFlux % mat (i,j) = ( waterInSnow % mat (i,j) - &
                                            swhc % mat (i,j) * swe % mat (i,j) ) / dtSnow
            !set maximum 
            waterInSnow % mat (i,j) = swhc % mat (i,j) * swe % mat (i,j)
            
        END IF
        
        
    END DO
END DO

!write point SWE
IF ( time == timePointExport ) THEN
   CALL SnowPointExport ( time )
END IF

RETURN
END SUBROUTINE SnowUpdate